
## Playing Politics with Environmental Protection: The Political Economy of Designating Protected Areas
## Jorge Mangonnet, Jacob Kopas, and Johannes Urpelainen
## August 2021

## File: functions for conducting regression analyses and producing tables and plots  
## NOTE: models will produce a list of regression objects that can then be fed into table functions


### REGRESSION MODELS ----------------------------------------------------------

## Main models ----
main.model <- function(formulas, code, data){
  for (i in 1:length(formulas)){
    model <-  felm(formula = formulas[[i]],
                   data = data)
    if (i == 1){
      model.list <- list(model)
    } else {
      model.list[[i]] <- model
    }
  }
  eval(parse(text = paste0(code, "<<- model.list")))
}

## Subset models (for comparing different subsets of the data) ----
sub.model <- function(formulas, code, datasets){
  for (j in 1:length(datasets)){
    data <- datasets[[j]]
    for (i in 1:length(formulas)){
      model <-  felm(formula = formulas[[i]],
                     data = data)
      if (i == 1 & j == 1){
        model.list <- list(model)
      } else {
        model.list[[length(model.list) + 1]] <- model
      }
    }
  }
  eval(parse(text = paste0(code, "<<- model.list")))
}



### REGRESSION TABLES ----------------------------------------------------------

omit.list <- NULL
omit.stats <- c("rsq", "ser", "f")
cov.labs <- NULL
title <- "Protected Areas"
order <- NULL
col.labels <- NULL
col.sep <- NULL

out.tab <- function(model.list, modeltype, ugrid){
  mpr <- unlist(lapply(X = model.list, function(x) length(unique(x$fe$m_pair))))
  spr <- unlist(lapply(X = model.list, function(x) length(unique(x$fe$s_pair))))
  if (modeltype == 0){
    add.lines <- list(c("Muni. Pair FE", rep("Yes", 6)),
                      c("Grid FE", "-", rep("Yes", 2), "-", rep("Yes", 2)),
                      c("State-Year FE", rep("-", 2), "Yes", rep("-", 2), "Yes"),
                      #c("Cluster SE", rep(c("Yes"), 6)),
                      c("Muni. Pairs", mpr),
                      c("Unique Grids", prettyNum(ugrid, ",")))
  } else if (modeltype == 1){
    add.lines <- list(c("Muni. Pair FE", rep("Yes", 3)),
                      c("Grid FE", "-", rep("Yes", 2)),
                      c("State-Year FE", rep("-", 2), "Yes"),
                      #c("Cluster SE", rep(c("Yes"), 3)),
                      c("Muni. Pairs", mpr),
                      c("Unique Grids", prettyNum(ugrid, ",")))
  } else if (modeltype == 2){
    add.lines <- list(c("Muni. Pair FE", rep("Yes", 9)),
                      c("Grid FE", rep(c("-", rep("Yes", 2)), 3)),
                      c("State-Year FE", rep(c(rep("-", 2), "Yes"), 3)),
                      #c("Cluster SE", rep(c("Yes"), 9)),
                      c("Muni. Pairs", mpr),
                      c("Unique Grids", prettyNum(ugrid, ",")))
  } else if (modeltype == 3){
    add.lines <- list(c("Muni. Pair FE", rep("Yes", 3)),
                      c("Grid FE",  rep("-", 3)),
                      c("State-Year FE", rep("-", 3)),
                      #c("Cluster SE", rep(c("Yes"), 3)),
                      c("Muni. Pairs", mpr),
                      c("Unique Grids", prettyNum(ugrid, ",")))
  } else if (modeltype == 4){
    add.lines <- list(c("Muni. Pair FE", rep("Yes", 6)),
                      c("Grid FE",  rep("-", 6)),
                      c("State-Year FE", rep("-", 6)),
                      #c("Cluster SE", rep(c("Yes"), 6)),
                      c("Muni. Pairs", mpr),
                      c("Unique Grids", prettyNum(ugrid, ",")))
  } else if (modeltype == 5){
    add.lines <- list(c("Muni. Pair FE", rep("Yes", 6)),
                      c("Muni. FE",  rep(c("-", "Yes", "Yes"), 2)),
                      c("State-Year FE", rep(c("-", "-", "Yes"), 2)),
                      #c("Cluster SE", rep(c("Yes"), 6)),
                      c("Muni. Pairs", mpr),
                      c("Unique Border Seg.", prettyNum(ugrid, ",")))
  } else if (modeltype == 6){
    add.lines <- list(c("Muni. Pair FE", rep("Yes", 4)),
                      c("Grid FE", "-","-", "Yes", "Yes"),
                      c("State-Year FE", "-","-", "-", "Yes"),
                      #c("Cluster SE", rep(c("Yes"), 4)),
                      c("Muni. Pairs", mpr),
                      c("Unique Border Seg.", prettyNum(ugrid, ",")))
  } else if (modeltype == 7){
    add.lines <- list(c("Muni. Pair FE", rep("Yes", 3)),
                      c("Grid FE", rep("Yes", 3)),
                      c("State-Year FE", rep("Yes", 3)))
                      #c("Cluster SE", rep("Yes", 3)))
  } else if (modeltype == 8){
    add.lines <- list(c("Muni. Pair FE", rep("-", 3)),
                      c("Muni FE", "-", rep("Yes", 2)),
                      c("State-Year FE", rep("-", 2), "Yes"))
                      #c("Cluster SE", rep("Yes", 3)))
  } else if (modeltype == 9){
    add.lines <- list(c("Muni. Pair FE", rep("-", 6)),
                      c("Muni FE", "-", rep("Yes", 2), "-", rep("Yes", 2)),
                      c("State-Year FE", rep("-", 2), "Yes", rep("-", 2), "Yes"))
                      #c("Cluster SE", rep("Yes", 6)))
  } else if (modeltype == 10){
    add.lines <- list(c("State Pair FE", rep("Yes", 3)),
                      c("Grid FE", "-", rep("Yes", 2)),
                      c("Year FE", rep("-", 2), "Yes"),
                      #c("Cluster SE", rep("Yes", 3)), 
                      c("State Pairs SE", spr), 
                      c("Unique Grids", prettyNum(ugrid, ",")))
  }  else if (modeltype == 11){
    add.lines <- list(c("Muni. Pair FE", rep("Yes", 4)),
                      c("State-Year FE", rep("Yes", 4)),
                      #c("Cluster SE", rep("Yes", 4)),
                      c("Muni. Pairs", mpr),
                      c("Unique Grids", prettyNum(ugrid, ",")))
  } else if (modeltype == 12){
    add.lines <- list(c("Covariates", rep("Yes", 4)),
                      c("Muni. Pair FE", rep("Yes", 4)),
                      c("Grid FE",  rep("-", 4)),
                      c("State-Year FE", "-", "Yes", "-", "Yes"),
                      #c("Cluster SE", rep(c("Yes"), 4)),
                      c("Muni. Pairs", mpr),
                      c("Unique Grids", prettyNum(ugrid, ",")))
  } 
    
  stargazer(model.list,
            header = F,
            float = F,
            no.space = T,
            covariate.labels = cov.labs,
            omit.stat = omit.stats,
            omit = omit.list,
            add.lines = add.lines,
            title = title,
            order = order,
            column.labels = col.labels,
            dep.var.labels = dep.labels,
            dep.var.labels.include = FALSE,
            star.cutoffs = c(0.10, 0.05, 0.01, 0.001),
            star.char = c("+", "*", "**", "***"),
            notes = c("$^{+}$p<0.1; $^{*}$p<0.05; $^{**}$p<0.01; $^{***}$p<0.001"),
            notes.append = FALSE,
            column.separate = col.sep)
}


### PLOTS-----------------------------------------------------------------------


## P-value plot ----
## NOTE: for covariate balance test (Figure A3)
plot.pval <- function(results, title=NULL, legend,legendx=0.15,legendy=2.2, textsize=0.9, parcex=0.8, at1=-0.35, at2=-0.15, at3=-0.9,xlim1=-0.85) {
  
  # set values of different parameters
  xlim = c(xlim1,1); pchset = c(21,24,22,23); pchcolset = c("blue","yellow","red","darkgreen")
  
  # set margins and letter size
  par(cex=parcex, mai = c(0.5, 0.35, 1.1, 0.35))
  
  # set number of rows 
  ny = nrow(results)
  
  # create the empty figure
  if(!is.null(title))  plot(x=NULL,axes=F, xlim=xlim, ylim=c(1,ny),xlab="",ylab="", main=title)
  if(is.null(title))   plot(x=NULL,axes=F, xlim=xlim, ylim=c(1,ny),xlab="",ylab="")
  
  # add the 0, 0.05 and 0.1 vertical lines
  abline(v=c(0,0.05,0.1),lty=c(1,4,4), lwd=c(1,2,2))
  axis(side=1,at=c(0,0.05,0.1,1),tick=TRUE, las=2, cex.axis=0.7)
  
  # add labels on top of the three areas of the graph
  axis(side=3,at=at1,labels="Mean\nTreated",tick=FALSE, padj=0.5,cex.axis=textsize)
  axis(side=3,at=at2,labels="Mean\nControl",tick=FALSE, padj=0.5,cex.axis=textsize)
  axis(side=3,at=0.5,labels="P-values",tick=FALSE, padj=0.5,cex.axis=textsize)
  
  # Fill the figure with the information which is inside the 'results' matrix
  # First, add the p-values as points
  for(i in 4:ncol(results)) points(results[,i],ny:1, pch = pchset[i-4+1], col = pchcolset[i-4+1], bg = pchcolset[i-4+1])
  
  # Second, add each variable name and the means for treated and control
  for(i in 1:ny) {
    text(at3,ny-i+1,results[i,1],adj = 0,cex=textsize) # variable name
    text(at1,ny-i+1,results[i,2], cex=textsize)        # treatment mean
    text(at2,ny-i+1,results[i,3], cex=textsize)        # control mean
  }
  
  # Add dotted horizontal lines every two variables to make it prettier
  for(i in seq(2,by=2,length.out=floor((ny-1)/2))) abline(h = i+0.5, lty = 3)
  
  # Add legend
  if(legend) legend(x=legendx, y=legendy, c("Diff-means", "Diff-medians", "D-stat", "Rank-sum"), pch=pchset, pt.bg = pchcolset, cex=0.8)
}


## Marginal effect plots (interactions) ----
# NOTE: for Figures 3, A7, A8, A9, A10, A12, and A13
marg.effect <- function(models, Z, order, labs){
  for (i in 1:length(models)){
    
    # Interaction terms and Coef estimates
    B1 <- models[[i]]$coefficients[order[1]]
    B3 <- models[[i]]$coefficients[order[2]]
    varB1 <- models[[i]]$clustervcv[order[1], order[1]]
    varB3 <- models[[i]]$clustervcv[order[2], order[2]]
    cvarB1B3 <- models[[i]]$clustervcv[order[1], order[2]]
    variable <- labs[i]
    marg <- B1 + (B3 * Z)
    
    # Variance
    m_se <- sqrt(varB1 + (Z^2 * varB3) + (2 * Z * cvarB1B3))
    
    m_lower <- marg - (1.96 * m_se)
    m_upper <- marg + (1.96 * m_se)
    out <- data.frame(x = Z, y = marg, m_lower = m_lower, m_upper = m_upper,
                      variable = variable)
    if (i == 1){
      m_data <- out
    } else {
      m_data <- rbind(m_data, out)
    }
  }
  
  # Output
  mar_eff <- ggplot(m_data, aes(x, y)) + geom_ribbon(aes(ymin=m_lower, ymax=m_upper), fill = "gray75") + 
    geom_line() + geom_line(y = 0, linetype = "dashed", color = "grey35") + geom_point(shape=".") +
    xlab("") + ylab("Marginal Treatment Effect") 
  
  return(mar_eff)
  
}
